Image analysis with EBImage
Image analysis with EBImage
EBImage is a Bioconductor package that offers tools to process and analyze images.
Installation
To install EBImage, start R (version “3.5”) and enter:
#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("EBImage", version = "3.8")Loading
To load EBImage, use the following command:
library("EBImage")Reading
Images are loaded using the function readImage(), which takes as input a character vector of file names or URLs. The loaded images are stored into an Image object.
#files=list.files(system.file("images",package="EBImage"))
#sampleColor=readImage(system.file("images","sample-color.png",package="EBImage"))
#sample=readImage(system.file("images","sample.png",package="EBImage"))
#shapes=readImage(system.file("images","shapes.png",package="EBImage"))
vesuvio=readImage("vesuvio.png")
cells=readImage(system.file("images","cells.tif",package="EBImage"))
nuclei=readImage(system.file("images","nuclei.tif",package="EBImage"))
greyVesuvio=channel(vesuvio,'grey')Grayscale images
str(greyVesuvio)## Formal class 'Image' [package "EBImage"] with 2 slots
## ..@ .Data : num [1:851, 1:315] 0.639 0.642 0.648 0.648 0.641 ...
## ..@ colormode: int 0print(greyVesuvio)## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 851 315
## frames.total : 1
## frames.render: 1
##
## imageData(object)[1:5,1:6]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.6392157 0.6339869 0.6418301 0.6483660 0.6418301 0.6392157
## [2,] 0.6418301 0.6418301 0.6431373 0.6457516 0.6379085 0.6379085
## [3,] 0.6483660 0.6522876 0.6418301 0.6392157 0.6496732 0.6562092
## [4,] 0.6483660 0.6418301 0.6392157 0.6483660 0.6509804 0.6457516
## [5,] 0.6405229 0.6339869 0.6444444 0.6483660 0.6418301 0.6339869dim(greyVesuvio)## [1] 851 315imageData(greyVesuvio)[10:14,4:8]## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.6483660 0.6444444 0.6366013 0.6444444 0.6483660
## [2,] 0.6444444 0.6366013 0.6405229 0.6444444 0.6405229
## [3,] 0.6444444 0.6405229 0.6483660 0.6522876 0.6444444
## [4,] 0.6483660 0.6483660 0.6444444 0.6562092 0.6483660
## [5,] 0.6483660 0.6444444 0.6444444 0.6405229 0.6366013Color images
print(vesuvio,short=TRUE)## Image
## colorMode : Color
## storage.mode : double
## dim : 851 315 3
## frames.total : 3
## frames.render: 1Image stacks
print(nuclei,short=TRUE)## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 510 510 4
## frames.total : 4
## frames.render: 4An image histogram shows the distribution of pixel intensities.
hist(greyVesuvio)
range(greyVesuvio)## [1] 0 1Displaying
The function display() allows to display images using an interactive JavaScript viewer or R’s built-in graphics capabilities.
display(greyVesuvio)
display(vesuvio)
display(nuclei,all=TRUE)
Writing
writeImage(greyVesuvio,'vesuvio.jpeg',quality = 85)
file.size('vesuvio.jpeg')## [1] 54781Image processing
The operations that can be applied to an image can be classified into three categories: point, local and global.
Point operations: output image pixel value at a specific position depends only on the input image pixel value at the same position. Examples: brigthness and contrast adjustment, gamma correction, thresholding.
Local operations: output image pixel value at a specific position depends on the pixel values in the neighborhood of input image pixel value at the same position. Example: convolution.
Global operations: output image pixel value at a specific position depends on all pixel values in the input image.
Brightness adjustment
\(y=x+b\)
display(
combine(
vesuvio+0.2,
vesuvio,
vesuvio-0.2
), all=TRUE
)
Contrast adjustment
\(y=x*c\)
display(
combine(
vesuvio*0.2,
vesuvio,
vesuvio*2
), all=TRUE
)
Gamma correction
\(y=x^{\gamma}\)
display(
combine(
vesuvio^0.7,
vesuvio,
vesuvio^2
), all=TRUE
)
Thresholding
Thresholding is a segmentation method that converts a gray-scale image into a binary image. During this operation all pixels whose intesity values are above a threshold are set to a foreground value, while all the remaining pixels are set to a background value.
\(y = x \le t\) or \(y = x \ge t\)
otsuGreyVesuvio = greyVesuvio < otsu(greyVesuvio)
print(otsuGreyVesuvio)## Image
## colorMode : Grayscale
## storage.mode : logical
## dim : 851 315
## frames.total : 1
## frames.render: 1
##
## imageData(object)[1:5,1:6]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] FALSE FALSE FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE FALSE FALSE
## [4,] FALSE FALSE FALSE FALSE FALSE FALSE
## [5,] FALSE FALSE FALSE FALSE FALSE FALSEdisplay(
combine(
greyVesuvio,
otsuGreyVesuvio
),all=TRUE
)
Morphological operations
Morphological operations are usually applied on binary images.
Erosion: the value of the output pixel is the minimum value of all the pixels in the input pixel’s neighborhood.
Dilation: the value of the output pixel is the maximum value of all pixels in the input pixel’s neighborhood.
Opening: is an erosion followed by a dilation.
Closing: is a dilation followed by an erosion.
display(
combine(
otsuGreyVesuvio,
erode(otsuGreyVesuvio),
dilate(otsuGreyVesuvio),
opening(otsuGreyVesuvio),
closing(otsuGreyVesuvio)
), all=TRUE
)
fillHull() allows to fill holes in objects.
display(
combine(
closing(otsuGreyVesuvio),
fillHull(closing(otsuGreyVesuvio))
), all=TRUE
)
Convolution
Convolution is a mathematical operation on two functions to produce a third function. Convolution can be written in the symbolic form as follows:
\(g(x,y)\) is the filtered image, \(f(x,y)\) is the original image, \(w\) is the kernel.
| \(I_{1}\) | \(I_{2}\) | \(I_{3}\) | \(I_{4}\) | \(I_{5}\) | \(I_{6}\) | \(I_{7}\) | \(I_{8}\) | \(I_{9}\) |
| \(I_{10}\) | \(I_{11}\) | \(I_{12}\) | \(I_{13}\) | \(I_{14}\) | \(I_{15}\) | \(I_{16}\) | \(I_{17}\) | \(I_{18}\) |
| \(I_{19}\) | \(I_{20}\) | \(I_{21}\) | \(I_{22}\) | \(I_{23}\) | \(I_{24}\) | \(I_{25}\) | \(I_{26}\) | \(I_{27}\) |
| \(I_{28}\) | \(I_{29}\) | \(I_{30}\) | \(I_{31}\) | \(I_{32}\) | \(I_{33}\) | \(I_{34}\) | \(I_{35}\) | \(I_{36}\) |
| \(I_{37}\) | \(I_{38}\) | \(I_{39}\) | \(I_{40}\) | \(I_{41}\) | \(I_{42}\) | \(I_{43}\) | \(I_{44}\) | \(I_{45}\) |
| \(I_{46}\) | \(I_{47}\) | \(I_{48}\) | \(I_{49}\) | \(I_{50}\) | \(I_{51}\) | \(I_{52}\) | \(I_{53}\) | \(I_{54}\) |
| \(I_{55}\) | \(I_{56}\) | \(I_{57}\) | \(I_{58}\) | \(I_{59}\) | \(I_{60}\) | \(I_{61}\) | \(I_{62}\) | \(I_{63}\) |
| \(I_{64}\) | \(I_{65}\) | \(I_{66}\) | \(I_{67}\) | \(I_{68}\) | \(I_{69}\) | \(I_{70}\) | \(I_{71}\) | \(I_{72}\) |
| \(I_{73}\) | \(I_{74}\) | \(I_{75}\) | \(I_{76}\) | \(I_{77}\) | \(I_{78}\) | \(I_{79}\) | \(I_{80}\) | \(I_{81}\) |
| \(K_{1}\) | \(K_{2}\) | \(K_{3}\) |
| \(K_{4}\) | \(K_{5}\) | \(K_{6}\) |
| \(K_{7}\) | \(K_{8}\) | \(K_{9}\) |
\(O_{12}=I_{2}K_{1}+I_{3}K_{2}+I_{4}K_{3}+I_{11}K_{4}+I_{12}K_{5}+I_{13}K_{6}+I_{20}K_{7}+I_{21}K_{8}+I_{22}K_{9}\)
Build kernel
Low-pass filter
A low-pass filter is used to blur images, remove artefacts. An example of low-pass filter is the gaussian filter. A filter can be generated using makeBrush().
gaussianKern=makeBrush(size = 5, shape = 'gaussian', sigma = 3)| 0.0317564 | 0.0375158 | 0.0396589 | 0.0375158 | 0.0317564 |
| 0.0375158 | 0.0443196 | 0.0468515 | 0.0443196 | 0.0375158 |
| 0.0396589 | 0.0468515 | 0.0495280 | 0.0468515 | 0.0396589 |
| 0.0375158 | 0.0443196 | 0.0468515 | 0.0443196 | 0.0375158 |
| 0.0317564 | 0.0375158 | 0.0396589 | 0.0375158 | 0.0317564 |
High-pass filter
A high-pass filter is used to detect edges, sharpen images. An example of high-pass filter is the laplacian filter.
values=c(-1,0,-1,0,4,0,-1,0,-1)
laplacianKern=matrix(values,nrow = 3,ncol = 3, byrow = TRUE)| -1 | 0 | -1 |
| 0 | 4 | 0 |
| -1 | 0 | -1 |
Laplacian of gaussian filter
values=c(0,0,-1,0,0,0,-1,-2,-1,0,-1,-2,16,-2,-1,0,-1,-2,-1,0,0,0,-1,0,0)
lapgauss5kern=matrix(values,nrow = 5,ncol = 5, byrow = TRUE)| 0 | 0 | -1 | 0 | 0 |
| 0 | -1 | -2 | -1 | 0 |
| -1 | -2 | 16 | -2 | -1 |
| 0 | -1 | -2 | -1 | 0 |
| 0 | 0 | -1 | 0 | 0 |
Filtering
filter2() allows to filter an image using the 2D FFT (2-Dimensional Fast Fourier Transform) convolution product.
GaussImage=filter2(vesuvio,gaussianKern)
LaplacianImage=filter2(vesuvio,laplacianKern)
LoGImage=filter2(vesuvio,lapgauss5kern)display (
combine(
vesuvio,
GaussImage,
LaplacianImage,
LoGImage
), all=TRUE
)
display (
combine(
channel(LaplacianImage,'grey'),
channel(LoGImage,'grey')
), all=TRUE
)
Example: application in cell biology
“cells.tif” and “nuclei.tif” are fluorescent microscopy images from two channels of HeLa cells, distributed with EBImage.
cells=readImage(system.file("images","cells.tif",package="EBImage"))
nuclei=readImage(system.file("images","nuclei.tif",package="EBImage"))
display(cells,all=TRUE)
display(nuclei,all=TRUE)
rgb = rgbImage(green=1.5*cells,blue=nuclei)
display(rgb,all=TRUE)
Nuclei segmentation
nmaskt = thresh(nuclei, w = 15, h = 15, offset = 0.05)
display(nmaskt,all=TRUE)
nmaskf = fillHull(opening(nmaskt, makeBrush(5, shape='disc')) )
display(nmaskf,all=TRUE)
display( combine(getFrame(nmaskt, 3), getFrame(nmaskf, 3)), all=TRUE )
dmap = distmap(nmaskf)
display(dmap,all=TRUE)
range(dmap)## [1] 0.00000 14.31782display(normalize(dmap), frame = 3)
nmask = watershed(dmap, tolerance = 2)
display(nmask,all=TRUE)
display( combine(
toRGB( getFrame(nuclei, 3) ),
colorLabels( getFrame(nmask, 3) )
), all=TRUE )
Cytoplasm segmentation
cmaskt = closing( gblur(cells, 1) > 0.105, makeBrush(5, shape='disc') )cmask = propagate(cells, seeds=nmask, mask=cmaskt, lambda = 0.001)display( combine(
toRGB( getFrame(cmaskt, 3) ),
colorLabels( getFrame(cmask, 3) )
), all=TRUE )
display( paintObjects(nmask,
paintObjects(cmask, rgb, col = "magenta", thick = TRUE),
col = "yellow", thick = TRUE), all = TRUE)
st = stackObjects(cmask, rgb)
display(st, all = TRUE)
features=computeFeatures.shape(cmask[,,1], cells[,,1])head(features)## s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
## 1 4216 408 43.52530 14.254634 16.13319 73.48543
## 2 3079 263 34.60167 9.061043 19.14618 53.49633
## 3 3088 211 30.99363 4.630083 16.81697 38.91545
## 4 2377 209 29.37053 8.030861 16.27813 44.85961
## 5 2643 277 34.75310 13.817002 15.03601 62.17141
## 6 2095 181 26.82256 6.948619 16.26868 40.91585
